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Abstract The case of the planar circular restricted three-body problem where one of the two 
primaries is an oblate spheroid is investigated. We conduct a thorough numerical analysis on 
the phase space mixing by classifying initial conditions of orbits and distinguishing between 
three types of motion: (i) bounded, (ii) escape and (iii) collisional. The presented outcomes 
reveal the high complexity of this dynamical system. Furthermore, our numerical analysis 
shows a strong dependence of the properties of the considered escape basins with the total 
orbital energy, with a remarkable presence of fractal basin boundaries along all the escape 
regimes. Interpreting the collisional motion as leaking in the phase space we related our 
results to both chaotic scattering and the theory of leaking Hamiltonian systems. We also 
determined the escape and collisional basins and computed the corresponding escape/crash 
times. The highly fractal basin boundaries observed are related with high sensitivity to initial 
conditions thus implying an uncertainty between escape solutions which evolve to different 
regions of the phase space. We hope our contribution to be useful for a further understanding 
of the escape and crash mechanism of orbits in this version of the restricted three-body 
problem. 

Keywords Restricted three-body problem; Escape dynamics; Escape basins; Eractal basin 
boundaries 


1 Introduction 


The issue of escape in Hamiltonian systems is a classical problem in nonlinear dynamics 
(e.g., Contopoulos 1 1990[l; [Contopoulos & Kaufmann] ( |1992| ); |Contopoulos et ah] ( [1993| ); 


[Schneider et al. 1 [2002 1 ). For energy levels above the escape energy the equipotential 
faces are open and exit channels emerge through which the particles can escape to infinity. 
The literature is replete with studies of such “open" Hamiltonian systems (e.g., [Barrio et| 
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ar|l |2009[l; Ernst & Peters] \20\A ) \ [Kandrup et al.| ( |1999| l; |Navarro & Henrard] ( |2001| l; |Zo-| 


tos| 1 2014a|b l). Nevertheless, the issue of escaping orbits in Hamiltonian systems is by far 


less explored than the closely related problem of chaotic scattering. In this situation, a test 
particle coming from infinity approaches and then scatters off a complex potential. This phe¬ 
nomenon is well investigated as well interpreted from the viewpoint of chaos theory (e.g., 
Bleher et al.|([7989ll;|Benet et al.|(|1996||1998|l;|Jung et al.|((T999|[7995)1 ;| Jung &~TeI|([T99T|l; 


Seoane et al. 


|2006[|; Seoane & Sanjuan 1 2007|l). At this point we should emphasize, that 


all the above-mentioned references on escapes and chaotic scattering in Hamiltonian sys¬ 
tem are exemplary rather than exhaustive, taking into account that a vast quantity of related 
literature exists. 

In open Hamiltonian systems an issue of paramount importance is the determination 
of the basins of escape, similar to basins of attraction in dissipative systems or even the 
Newton-Raphson fractal structures. An escape basin is defined as a local set of initial con¬ 
ditions of orbits for which the test particles escape through a certain exit in the equipotential 
surface for energies above the escape value. Basins of escape have been studied in many 
earlier papers (e.g., [Bleher et al.| (| 198 8[l; [Kennedy & Yorke|(|1991 1 ; Poon et al. l[1996| l]). 


The reader can find more details regarding basins of escape in (|Contopoulos[ 2002^, while 


the review ( |Zotos[ |2014b^ provides information about the escape properties of orbits in a 
multi-channel dynamical system of a two-dimensional perturbed harmonic oscillator. The 
boundaries of an escape basins may be fractal (e.g., [Aguirre et aT (|2009^; Bleher et al.j 


([1988[l) or even respect the more restrictive Wada property (e.g., Aguirre et al. (|2001 1 ), 


the case where three or more escape channels coexist in the equipotential surface. Escaping 
orbits in the classical Restricted Three-Body Problem (RTBP) is another typical example 
(e.g., [Nagler[[2004|[2005[[de Assis & Terra[[2014[ l. 

The classical RTBP assumes that the masses of the two primaries are spherically sym¬ 
metrical in homogeneous layers however, it is found that several celestial bodies, such as 


Saturn and Jupiter are sufficiently oblate (|Beatty et a^ 

1999|l. In addition, the minor planets 

(e.g., Ceres) and meteoroids have irregular shapes (Mi 

lis et al.|[1987[[Norton & Chitwood[ 


[2008[ l. The oblateness or triaxiality of a celestial body can produce perturbation deviations 
from the two-body motion. The most striking example of perturbations arising from oblate¬ 
ness in the solar system is the orbit of the fifth satellite of Jupiter, Amalthea. This planet is 
so oblate and the satellite’s orbit is so small that its line of apsides advances about 900° in a 


year (e.g., Moulton ([1914 1 ). The study of oblateness includes the series of works of Beevi & 

Sharma[(|2012|);[Markellos et al.[(|2002|l; 

Kalantonis et al.[([2005[[20061[20081l; 

Kalvouridis & 

Gousidou-Koutita|(|2012J;[Perdiou et al 

(|2012^;[Sharma & Subba Rao|(11979 

198611; Subba 

Rao & Sharma ([1988 1997|l; Sharma 

1981 [1987 1989 1990^ by considering the more 


massive primary as an oblate spheroid with its equatorial plane co-incident with the plane 
of motion of the primaries. 

It is well known that all the planets of the Solar System are not perfect spheres but 
spheroidals. Therefore when someone desires to model a three-body system it is essential 
to properly modify the RTBP in order to take into account the exact shape of the planets 
(primary bodies). In [Oberti & Vienn^ ( |2003[ l the authors compare observational data from 
the systems Saturn-Tethys-satellite and Saturn-Dione-satellite and they conclude that the 
corresponding theoretical data are much more accurate when the oblateness of Saturn is 
taken into consideration. In the same vein, in the work of [Stuchi et al.[ ( |^08[ l regarding the 
dynamics of a spacecraft in the Neptune-Triton problem the oblateness coefficient is also 
included. In the present paper we continue the work initiated in [Nagler] ( |2004[ ) and [Nagler| 
( [2005 [ I (hereafter Paper I and II) following the same numerical techniques. Our aim is to 
numerically investigate the properties of motion in the RTBP with oblateness. As far as we 
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Fig. 1 Schematic picture of the planar circular restricted three-body problem for the Copenhagen case. 


know, this is the first detailed and systematic numerical analysis on the phase space mixing 
of bounded motion, escape and crash in the RTBP with oblateness and this is exactly the 
novelty and the contribution of the current work. 

The paper is organized as follows: In Section|^we introduce the considered dynamical 
model and we present its properties along with some necessary details. All the computa¬ 
tional methods we used in order to determine the character of orbits are described in Section 
In the following Section, we conduct a thorough numerical investigation revealing the 
overall orbital structure (bounded regions and basins of escape/crash) of system and how 
it is affected by the total orbital energy considering three cases regarding the value of the 
oblateness coefficient. In Section]^ a general overview is provided showing in more detail 
the influence of the energy as well as the oblateness parameter. Our paper ends with Section 
1^ where the discussion and the conclusions of this research are given. 


2 Details of the dynamical model 

Let us briefly recall the basic properties and some aspects of the planar restricted threeaA§- 
body problem ( |SzebeheIyl|1967[ l. The two primaries move on circular orbits with the same 
Kepler frequency around their common center of gravity, which is assumed to be fixed at 
the origin of the coordinates. The third body (test particle with mass much smaller than 
the masses of the primaries) moves in the same plane under the gravitational field of the 
two primaries (see Fig.[^. The non-dimensional masses of the two primaries are 1 — jU and 
jU, where /i = m 2 /{m\ + m 2 ) is the mass ratio. We consider the Copenhagen case where 
M = 1/2. 

We choose as a reference frame a rotating coordinate system where the origin is at (0,0), 
while the centers Ci and C 2 of the two primaries are located at (—/i,0) and (1 —/i,0), 
respectively. The total time-independent gravitational potential is ([Sharma & Subba Rao| 
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ri n 2r{ 2 ' 

where 

n = \/ {x + ^lf+y^, r2 = \J {x + ^-\f +y^, w^ = l + ^, 


( 1 ) 

(2) 


are the distances to the respective primaries and the angular velocity (n), while Ai is the 
oblateness coefficient which is defined as 


Ai 


{RE)^ - {RPf 

5^2 


(3) 


where RE and RP are the equatorial and polar radius, respectively of the oblate primary, 
while R is the distance between the centers of the two primaries. We consider values of the 
oblateness coefficient is in the interval [0,0.1] (see e.g., (Kalantonis et al. 2006| Perdios 


|& Kalantonis ] 2006 1 ), while we study the effect of oblateness up to the linear coefficient 
J 2 only ( Abouelmagd[ |2012] l . We must point that the values of the oblateness in the Solar 
system are relatively low (e.g. Table 1 in |Sharma & Subba Rao| l |T976[ l). 

The scaled equations of motion describing the motion of the test body in the corotating 
frame read ( jSharma & Subba Rao||1976| l 


.. . dV{x,y) . dV{x,y) 

x = 2ny -r-, y=—2nx --. 

ox oy 

The dynamical system 0 admits the well know Jacobi integral 


J{x,y,x,y} = - (x^+f) +V{x,y)=E, 


(4) 

(5) 


where x and y are the momenta per unit mass, conjugate to x and y, respectively, while E is 
the numerical value of the energy which is conserved and defines a three-dimensional invari¬ 
ant manifold in the total four-dimensional phase space. Thus, an orbit with a given value of 
it’s energy integral is restricted in its motion to regions in which E <V {x,y), while all other 
regions are forbidden to the test body. It is widely believed that J is the only independent 
integral of motion for the PCRTBP system ( |Poincare]|1993| l. The energy value E is related 
with the Jacobi constant by C = —2E. 

The dynamical system has five equilibria known as Lagrangian points l [SzebeheI^|1967^ 
at which 


dy{x,y) ^ dv{x,y) 

dx dy 


( 6 ) 


The isolines contours of constant potential, the position of the five Lagrangian points Z,, , i = 
1,5, as well as the centers of the two primaries are shown inFig.j^whereAi = 0.01. Three of 
them, El, L 2 , and L 3 , are collinear points located in the jc-axis. We note here that the isolines 
contours are symmetrical only with respect to the y = 0 axis, where the symmetry to the 
X = 0 axis is lost due to the oblateness. The central stationary point L\ is a local minimum 
of the potential V{x,y). The stationary points L 2 and L 3 are saddle point. Let E 2 located 
at X < 0, while L 3 be at x > 0. The points and L 5 on the other hand, are local maxima 
of the gravitational potential, enclosed by the banana-shaped isolines. The projection of 
the four-dimensional phase space onto the physical (or position) space (x,y) is called the 
Hill’s regions and is divided into three domains shown in Fig. with different colors; (i) the 
interior region (green) for x(L 2 ) < ^ (ii) the exterior region (yellow) for x < x{L2) 
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Fig. 2 The isolines contours of the constant potential, the location of the centers of the two primaries (blue) 
and the position of the five Lagrangian points (red), for A; = 0.01. The interior region is indicated in green, 
the exterior region is shown in yellow, while the forbidden regions of motion are marked with grey. 


and X > x{L-i)\ (iii) the forbidden regions (gray). The boundaries of these Hill’s regions are 
called Zero Velocity Curves (ZVCs) because they are the locus in the physical (x,y) space 
where the kinetic energy vanishes. Finally, the position of the Lagrangian points L 2 and L 3 
is a function of the oblateness coefficient Ai (e.g., |Singh & Lek^p012[ l). 

The values of the Jacobi integral at the five Lagrangian points L, are critical energy levels 
and are denoted as £,■ (Note that E\ = 0, while £4 = £5). The structure of the equipotential 
surfaces strongly depends on the value of the energy. In particular, there are four distinct 
cases 

- E < E 2 '. Both necks at L 2 and £3 are closed, so in the interior region we have only 
collisional and bounded motion. 

- E 2 <E <Ey. Only the neck around £3 is open which acts as escape channel. 

- E 2 < E < E 4 : The necks around both £2 and £3 are open and two symmetrical, with 
respect to the y = 0 axis, forbidden regions are present. 

- £ > £4: The banana-shaped forbidden regions disappear and therefore, motion over the 
entire physical (x,y) plane is possible. 

In Fig. [3 a-d) we present for Ai = 0.01 a characteristic equipotential surface for the four 
possible Hill’s region configurations. We observe in Fig.[^ the two openings (exit channels) 
at the Lagrangian points £2 and £3 through which the body can leak out. In fact, we may 
say that these two exits act as hoses connecting the interior region of the system where 
x(£ 2 ) <x < ^(£ 3 ) with the “outside world" of the exterior region. 


3 Computational methods and criteria 

The motion of the test third body is restricted to a three-dimensional surface £ = const, 
due to the existence of the Jacobi integral. With polar coordinates {r,(j)) in the center of 
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Fig. 3 Four possible Hill’s region configurations for the PCRTBP system when A; = 0.01. The white domains 
correspond to the Hill’s region, gray shaded domains indicate the forbidden regions, while the thick black 
lines depict the Zero Velocity Curves (ZVCs). The red dots pinpoint the position of the Lagrangian points, 
while the positions of the centers of the two primaries are indicated by blue dots, (a-upper left): E = —1.76; 
(b-upper right): E = —1.742; (c-lower left): E = —1.73; (d-lower right): E = —1.30. 


the mass system of the corotating frame the condition f = 0 defines a two-dimensional 
surface of section, with two disjoint parts (j) < 0 and (j) > 0. Each of these two parts has a 
unique projection onto the configuration physical (jc,y) space. Our investigation takes place 
in both types of projection for a better understanding of the orbital dynamics. In order to 
explore the behavior of test particles in the Copenhagen model, we need to define samples 
of initial conditions of orbits whose properties will be identified. Fore this purpose, we define 
for several values of the total orbital energy E, dense uniform grids of 1024 x 1024 initial 
conditions regularly distributed on the (jc,y) plane inside the area allowed by the value of the 
energy. Following a typical approach, the orbits are launched with initial conditions inside a 
certain region, called scattering region, which in our case is a square grid with —2 < x,y < 2. 

In the PCRTBP system the configuration space extends to infinity thus making the iden¬ 
tification of the type of motion of the test body for specific initial conditions a rather de¬ 
manding task. There are three possible types of motion for the test body: (i) bounded motion 
around one of the primaries, or even around both; (ii) escape to infinity; (iii) crash into one 
of the primaries. Now we need to define appropriate numerical criteria for distinguishing 
between these three types of motion. The motion is considered as bounded if the test body 
stays confined for integration time fmax inside the system’s disk with radius and cen¬ 
ter coinciding with the center of mass origin at (0,0). Obviously, the higher the values of 
fmax and Rd the more plausible becomes the definition of bounded motion and in the limit 
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Fig. 4 Schematic picture of the three different types of motion. The motion is considered to be bounded if 
the test body stays confined for integration time fmax inside the system’s disk with radius Rj = 10, while the 
motion is unbounded and the numerical integration stops when the test body crosses the system’s disk with 
velocity pointing outwai'ds. Crash with one of the primaries occurs when the test body crosses the disk of 
radius R,„i = S ,„2 = 10^“* of one of the primaries. 


fmax °° the definition is the precise description of bounded motion in a finite disk of radius 
Rj. Consequently, the higher these two values, the longer the numerical integration of initial 
conditions of orbits lasts. In our calculations we choose tmax = 10^ and Rj = 10 (see Fig. 
1^. We decided to include a relatively high disk radius {Rd = 10) in order to be sure that 
the orbits will certainly escape from the system and not return back to the interior region. 
Furthermore, it should be emphasized that for low values of t^ax the fractal boundaries of 
stability islands corresponding to bounded motion become more smooth. Moreover, an orbit 
is identified as escaping and the numerical integration stops if the test body body intersects 
the system’s disk with velocity pointing outwards at a time tesc < tmax- Finally, a crash with 
one of the primaries occurs if the test body, assuming it is a point mass, crosses the disk 
with radius around the primary, where in our case we choose = 10^^. Here is should 
be noted that in generally it is assumed that the radius of a celestial body (e.g., a planet) is 
directly proportional to the cubic root of its mass. For the sake of simplicity of the numer- 

( |2006[ l). In papers I and II it was shown that the radii of the primaries influence the area of 
crash and escape basins. 

The vast majority of bounded motion corresponds to initial conditions of regular orbits. 
It therefore seems appropriate to further classify initial conditions of ordered orbits into 
regular families. For this task we use the symbolic orbit classification which was also used 
in Papers I and II. According to this method orbits are classified by taking into account 
their orientation with respect to the centers of the two primaries C\ and C 2 , as well as their 
rotation (clockwise or counterclockwise). In particular, the orbit classification is based on an 
automatic detection of x axis passages of the test body. Furthermore, two consecutive x axis 
passages define a half rotation with respect to the fixed centers of the two primary bodies. 


ical calculations we decided to fix the radii of the primaries to 10 ^ (see also 


Barrio et al. 
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Fig. 5 Characteristic orbit examples of the seven main types of regular orbits. 


In Fig. we present characteristic orbit examples of the seven main types of regular orbits. 
A more precise description of the types of orbits can be found in |Nagler1 ( [2002| ). 

As it was stated earlier, in our computations, we set 10^ time units as a maximum time 
of numerical integration. The vast majority of escaping orbits (regular and chaotic) however, 
need considerable less time to escape from the system (obviously, the numerical integration 
is effectively ended when an orbit moves outside the system’s disk and escapes). Neverthe¬ 
less, we decided to use such a vast integration time just to be sure that all orbits have enough 
time in order to escape. Remember, that there are the so called “sticky orbits" which behave 
as regular ones during long periods of time. Here we should clarify, that orbits which do 
not escape after a numerical integration of 10^ time units are considered as non-escaping or 
trapped. 

The equations of motion Q for the initial conditions of all orbits are forwarded in- 
tegrate d using a double precision Bulirsch-Stoer FORTRAN 77 algorithm (e.g., [Press et ak] 
1 1992 1 ) with a small time step of order of 10^^, which is sufficient enough for the desired 
accuracy of our computations. Here we should emphasize, that our previous numerical ex¬ 
perience suggests that the Bulirsch-Stoer integrator is both faster and more accurate than a 
double precision Runge-Kutta-Fehlberg algorithm of order 7 with Cash-Karp coefficients. 
Throughout all our computations, the Jacobian energy integral (Eq. l|^) was conserved bet¬ 
ter than one part in 10^*^, although for most orbits it was better than one part in 10^*^. For 
collisional orbits where the test body moves inside a region of radius 10^^ around one of 
the primaries the Lemaitre’s global regularization method is applied. 


4 Numerical results & Orbit classification 

The main objective of our investigation is to classify initial condition of orbits in the physical 
(jc,y) plane into three categories: (i) bounded orbits; (ii) escaping orbits and (iii) crashing 
orbits, distinguishing simultaneously regular orbits into different types. Furthermore, two 
additional properties of the orbits will be examined: (i) the time-scale of crash and (ii) the 
time-scale of the escapes (we shall also use the terms escape period or escape rates). In the 
present paper, we shall explore these dynamical quantities for various values of the total 
orbital energy, as well as for the oblateness coefficient Ai. In particular, three different cases 
are considered: (a) the primary body with center at Ci has a small value of Ai, (b) the 
primary body has an intermediate value of Ai and (c) the primary body is highly oblate. 
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In the following color-coded grids (or orbit type diagrams - OTDs) each pixel is assigned 
a color according to the orbit type. Thus the initial conditions of orbits are classified into 
bounded motion of a few types, unbounded escaping motion and collisional motion. In this 
special type of Poincare surface of section the phase space emerges as a close and compact 
mix of escape basins, crash basins and stability islands. 

Our numerical calculations indicate that apart from the escaping and crashing orbits 
there is also a considerable amount of non-escaping orbits. In general terms, the majority of 
non-escaping regions corresponds to initial conditions of regular orbits, where an adelphic 
integral of motion is present, restricting their accessible phase space and therefore hinders 
their escape. Additional numerical computations reveal that the bounded regular orbits are 
mainly loop 1:1 resonant orbits for which the adelphic integral applie^ while other types 
of secondary resonant orbits are also present. The n : m notation we use for the regular orbits 
is according to |Carpintero & Aguilar| ( [T998[ l and |Zotos & Carpintero| ( [20T^ , where the ratio 
of those integers corresponds to the ratio of the main frequencies of the orbit, where main 
frequency is the frequency of greatest amplitude in each coordinate. Main amplitudes, when 
having a rational ratio, define the resonances of an orbit. 


4.1 Case I: Results for a low value of oblateness 

Our exploration begins considering the case where the primary body 1 has a relatively low 
value of oblateness, that is the case of Ai = 0.001. In Fig. the OTD decompositions for 
both (j) <0 (left column) and ^ > 0 (right column) reveal the structure of the physical {x,y) 
space for three energy levels, where the several types of orbits are indicated with different 
colors. The color code is explained in the color bar at the bottom of the figure. The black 
solid lines in the two types of plots denote the Zero Velocity Curve, while the inaccessible 
forbidden regions are marked in gray. The color of a point represents the orbit type of a 
test body which has been launched with pericenter position at {x,y). The three energy levels 
belong to each of the three Hill’s regions configurations explained earlier in Fig.[^ The cases 
E < E 2 and E 2 < E < Et, give very similar results, so we included only the latter. When 
E = —1.73 we observe that in both cases the interior region is filled with initial conditions 
of orbits that either are regular or crash into one of the primaries, while initial conditions 
of escaping orbits are present only in the exterior region outside L 2 and L 3 . Regular motion 
dominates the interior region and two large stability islands are shown near the centers of the 
primaries. These stability islands contain quasi-periodic orbits which are symmetrical with 
respect to a reflection over the x axis and they move in clockwise sense, hence, retrograde in 
relation to the rotating system of coordinates. Moreover, the stability regions are surrounded 
by a mix of domains of crash orbits with respect to the first and the second primary body. 
On the other hand, the exterior region contains mostly initial conditions of escaping orbits 
however, in the (j) <0 plot we have to point out the existence of a stability ring containing 
initial conditions of regular orbits that circulate clockwise around both primaries. As the 
value of the energy increases and the area of forbidden region is reduced it is seen that the 
stability ring disappears, while in both the interior and exterior regions basins of escaping 
and collisional orbits are formed. We should note though, that the basins containing orbits 
that crash into primary 2 are much smaller, with respect to the crash basins of primary 1 , 
and they appear only as thin filaments. This phenomenon is justified by the high value of 

* The total angular momentum is an approximately conserved quantity (integral of motion), even for orbits 
in non spherical potentials. 
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Fig. 6 The orbital structure of the physical (jc,y) plane in a corotating frame of reference is given using orbit 
type diagrams (OTDs) for three energy levels and for both parts (j) <0 (left column) and (j> >0 (right column) 
of the surface of section r = 0, when Ai — 0.001. (Top row): E = — 1.73; (middle row): E — — 1 .60; (bottom 
row): E — —1.30. The vertical black dashed lines denote the centers of the two primaries, wile the vertical 
purple dashed lines indicate the position of the Lagrangian points L 2 and L 3 . The color bar contains the color 
code which relates the types of orbits presented in Fig.[^with different colors. 
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Fig. 7 Distribution of the escape and collisional time of the orbits on the physical (jr,y) space when = 
0.001 for the energy levels of Fig.(Top row): E — —1.73; (middle row): E — —1.60; (bottom row): E — 
— 1.30. The darker the color, the larger the escape/crash time. Initial conditions of bounded regular orbits are 
shown in white. 
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Fig. 8 Evolution of the percentages of all types of orbits on the physical {x,y) plane, for Aj = 0.001 when 
varying the energy parameter E for (a-left): the 0 < 0 part and (b-right): the 0 > 0 part of the configuration 
space. 


the oblateness of primary 1. For £ = — 1.30 > £4 the test body has full access to the entire 
physical {x,y) plane. Two regions of bounded motion are shown around the primaries, where 
for each stability island there is a parent periodic orbit at its center. Furthermore, due to the 
rotation of the primaries the crash basins wind out in spiral form in the outer regions of the 
(0 < 0 plot. Crash basins are also present in the immediate vicinity of the origin. At this point 
we would like to stress out that the area of crash basins is several orders of magnitude larger 
than the total size of the primary body’s disks. In the OTD for (0 > 0 the escape basin cover 
the vast majority of the configuration space, while the total size of stability regions is less 
than for (0 < 0. It is interesting to note that all regular orbits in all three energy levels were 
found to be retrograde thus traveling in clockwise sense. 

The following Fig. shows how the escape and crash times of orbits are distributed 
on the physical {x,y) space for the three energy levels discussed in Fig. Light reddish 
colors correspond to fast escaping/crashing orbits, dark blue/purple colors indicate large 
escape/crash rates, while white color denote stability islands of regular motion. Note that the 
scale on the color bar is logarithmic. Inspecting the spatial distribution of various different 
ranges of escape time, we are able to associate medium escape time with the stable manifold 
of a non-attracting chaotic invariant set, which is spread out throughout this region of the 
chaotic sea, while the largest escape time values on the other hand, are linked with sticky 
motion around the stability islands of the two primaries. It should be noted that the behaviour 
of the escape times is very similar to that observed in |de Assis & Ter^ ( |2014[ l. As for the 
collision time we see that orbits with initial conditions very close to the vicinity of the center 
of the oblate primary 1 collide with it almost immediately, within the first time step of the 
numerical integration. 

The evolution of the percentages of the different types of orbits on the physical (jc,y) 
space for both parts (0 < 0 (a) and (0 > 0 (b) when the energy varies is presented in Fig. 
[^a-b). The vertical black dashed lines indicate the three critical values of the Jacobi inte¬ 
gral (£ 2 ,£ 3 ,£ 4 ). It is seen that for Ai = 0.001, £2 and £3 are very close and therefore the 
corresponding lines are almost indistinguishable. One may observe that in both cases three 
types of orbits control the vast majority of the configuration space: (i) escaping orbits; (ii) 













Crash test for the restricted three body problem with oblateness 


13 


orbits which crash into oblate primary 1 and (iii) type 3b regular orbits that is orbits which 
circulate clockwise around both primaries, while all the rest types are almost unaffected by 
the energy shifting, having considerably low rates (less than 5%) throughout. For the (j) <0 
case we see that for low values of energy E < E 2 , escaping and type 3b regular orbits seem 
to share about 90% of the physical space, while initial conditions of orbits that crash into 
oblate primary 1 occupy only about 7% of the same plane. As the value of the energy in¬ 
creases however, the percentage of orbits that crash into the oblate primary increases and for 
Z? > —1.3 it dominates, while the rate of type 3b regular orbits on the other hand, reduces 
drastically and for £ > —1.6 it completely vanishes. Moreover, the percentage of escaping 
orbits increases linearly for E < —1.6, while for larger values of the energy this tendency is 
reversed. At the highest energy level studied (£ = —1), about 35% of the physical space is 
occupied by escaping orbits, while initial conditions of orbits that lead to collision with pri¬ 
mary 1 correspond to about 60% of the {x,y) plane. For the the <j) > 0 case, things are much 
more simpler since always the vast majority (more than 80%) of the configuration space is 
dominated by escaping orbits, while the rate of collisional orbits to primary 1 is significantly 
lower; it starts at about 20% at low values of energy {E < £ 2 ) and drops to about 5% for 

£ > £4. 


4.2 Case II: Results for an intermediate value of oblateness 

We continue our numerical quest considering the case where the oblateness coefficient of 
the primary body 1 has an intermediate value, that is Ai = 0 . 01 , and we follow the same 
numerical approach as previously. The structure of the physical (x,y) space for three energy 
levels is unveiled in Fig. [^through the OTD decompositions for both (p <0 (left column) 
and (j) > 0 (right column) parts of the configuration space. We observe that in general terms 
for both cases {(j) <0 and <j) > 0), things are quite similar to those discussed earlier in Fig. 

However there is one major difference regarding the stability regions. In particular, for 
the energy level £ = —1.742, that is for the Hill’s region £2 < £ < £ 3 , the stability island 
near the center of the oblate primary 1 is absent. At £ = —1.6 a tiny stability region emerges 
at the left and right part of the center for (^ < 0 and (j) > 0, respectively, while only at a 
relatively high energy level £ = —1.3 > £4 bounded retrograde motion around primary 1 is 
possible. Therefore, we may conclude that the oblateness coefficient influences significantly 
the regular orbits around the oblate primary 1. Our numerical computations strongly indicate 
that the more oblate is primary 1 the less bounded motion is observed around it, while at the 
same time most of the initial conditions of orbits launched relatively close to primary 1 lead 
to fast collision with it. Moreover, as we seen in Fig. the (p > 0 part of the configuration 
space is dominated by escaping orbits, while inside the interior region crash with the left 
oblate primary body becomes more and more likely as we proceed to higher energy levels. 
Around the spherical (non-oblate) primary body 2 on the other hand, the area occupied 
by initial conditions corresponding to bounded motion around primary 2 exhibits a minor 
decrease with increasing energy. In addition, the crash basin to primary 2 is much smaller 
with respect to the extended crash basin 1, and has a spiral shape around primary 2. The 
distribution of the escape and crash times of orbits on both parts of the physical space is 
shown in Fig.[T^ One may observe that the results are very similar to those presented earlier 
in Fig.jTj where we found that orbits with initial conditions inside the escape and crash basins 
have the smallest escape/crash rates, while on the other hand, the longest escape/crash times 
correspond to orbits with initial conditions in the fractal regions of the plots. It is interesting 
to note that orbits with initial conditions in the vicinity of the center of the oblate primary 
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Fig. 9 The orbital structure of the physical (jc,y) plane in a corotating frame of reference is given using orbit 
type diagrams (OTDs) for three energy levels and for both parts (j) <0 (left column) and (j> >0 (right column) 
of the surface of section r = 0, when = 0.01. (Top row): E = —1.742; (middle row): E — —1.60; (bottom 
row): E — —1.30. The vertical black dashed lines denote the centers of the two primaries, wile the vertical 
purple dashed lines indicate the position of the Lagrangian points L 2 and L 3 . The color bar contains the color 
code which relates the types of orbits presented in Fig.[^with different colors. 
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Fig. 10 Distribution of the escape and collisional time of the orbits on the physical (A',y) space when Aj = 
0.01 for the energy levels of Fig. (Top row): E — — \ .742; (middle row): E — — \ .60; (bottom row): E — 
— 1.30. The darker the color, the larger the escape/crash time. Initial conditions of bounded regular orbits are 
shown in white. 
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Fig. 11 Evolution of the percentages of all types of orbits on the physical (a:,)’) plane, for Ai = 0.01 when 
varying the energy parameter E for (a-left): the 0 < 0 part and (b-right): the 0 > 0 part of the configuration 
space. 


1 collide with it almost immediately, while this phenomenon is not observed for the case of 
primary body 2 which is not oblate. 

In Fig. m a-b) we demonstrate the evolution of the percentages of the different types of 
orbits on the physical (x,y) space for both parts (j) <0 (a.) and (^ > 0 (b) as a function of the 
total orbital energy. We observe that for Ai =0.01 the vertical black dashed lines indicating 
the critical values of the energy E 2 and £3 are no longer indistinguishable as it was in Fig. 
[^for Ai = 0.001. Once more it is seen that in both parts of the configuration space only 
three types of orbits (escaping, crash into primary 1 and regular type 3b) are influenced by 
the change on the value of the energy, while all the other types of orbits seem completely 
unaffected by the energy shifting holding very low rates (less than about 5% throughout the 
energy range). In Fig. [nji we see that for low value of the energy E < E 2 the percentages of 
escaping and type 3b regular orbits seem to coincide at about 45%, while for larger energies 
they follow completely different paths. In particular, the rate of escaping orbits initially 
increases but for £ > £3 it exhibits a constant and almost linear decrease, while the rate of 
type 3b ordered orbits displays a rapid reduction and for £ > —1.6 it completely vanishes. 
The percentage of orbits that crash into primary 1 on the other hand, increases drastically 
and for £ > — 1.5 this type of orbits is the most populated family. Moreover, at the highest 
energy level studied, that is £ = — 1 , crashing into primary 1 orbits cover about two thirds 
of the entire physical {x,y) space. The orbital structure of the (j) > 0 part of the configuration 
space is completely different. Indeed, in Fig. [Ujt one may observe that escaping orbits 
dominate the physical plane throughout the energy range. For low values of the energy the 
corresponding rate ia about 80%, while for £ > £4 it exceeds 95%. The evolution of the 
percentage of orbits that crash into oblate primary 1 follows an opposite pattern, starting 
from about 20% and drops to about 5% for £ > —1.2. We should point out that in this case, 
the rates of all the other types of orbits are always less than 2 %, while the change on the 
value of the energy only shuffles the orbital content among them. 
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4.3 Case III: Results for a high value of oblateness 


The last case under investigation involves the scenario where the oblateness coefficient of the 
primary body 1 has a relatively high value, that is Ai =0.1. Again, all the different aspects 
of the numerical approach remain exactly the same as in the two previously studied cases. 
The orbital structure of the physical {x,y) plane through the OTD decompositions for both 
(j) <0 (left column) and (j) >0 (right column) parts of the configuration space distinguishing 
between the main types of orbits for three energy levels is shown in Fig. 12 For the (j) <0 


part it is evident that when E = —1.90 < C 3 all the interior region is occupied either by 
collisional orbits or regular orbits, while the exterior region is divided into two domains; a 
circular domain containing initial conditions of escaping orbits and another domain covered 
by type 3b ordered orbits that circulate clockwise around both primaries. The collisional to 
primary 2 basin is appeared as thin layer outside the main stability island. As the value of the 
energy increases the extent of this stability island is reduced, while a considerable amount 
of initial conditions of orbits located inside the interior region lead to escape. At the same 
time, a large portion of the exterior region is occupied by extended crash basins spiralling 
around the center of the coordinates. The (^ > 0 part of the configuration space displays, once 
more, the same structure as in the previous two cases, where we found that the vast majority 
of the integrated initial conditions correspond to orbits which escape from the system. In 
the previous case where the primary body 1 had an intermediate value of oblateness, we 
deduced that the more oblate is one primary the less bounded motion is observed around 
it. This conclusion becomes more clear and strong here where we see for a high value of 
the oblateness (Ai = 0 . 1 ) there is no indication of regular motion around primary 1 for all 
tested energy levels and in both parts of the configuration space. Thus we may certainly 
conclude that the increase on the oblateness of the primary body has a detrimental effect 
on the stability islands near the same primary, leading most test bodies in collisional orbits 
with the oblate primary. In Fig. we depict the distribution of the escape/crash times of 
orbits for both parts of the physical space, where one can see similar outcomes with that 
presented in the two previous subsections. At this point, we would like to point out that the 
basins of escape can be easily distinguished in Fig.jT^ being the regions with intermediate 
greenish colors indicating fast escaping orbits. Indeed, our numerical calculations suggest 
that orbits with initial conditions inside these basins need no more than 10 time units to 
escape from the system. Furthermore, the crash basins are shown with reddish colors where 
the corresponding crash time is less than one time unit. 

Finally, Fig. [T4|a-b) shows how the percentages of all types of orbits on both parts of 
the physical space evolve when the total orbital energy varies in the interval E 6 [—2, — 1]. 
In this case, the four energy intervals defined by the three critical values of the energy are 
fully distinguishable and are the following: (i) E < E 2 ', (ii) E 2 <E <Ey, (iii) E^ < E < £4; 
(iv) E > Eli. Fig- [14^ shows that in the first energy interval escaping and type 3b regular 
orbits occupy about 40% and 50% of the physical plane, respectively, while the rest area 
corresponds to initial conditions of orbits that crash into primary 1. For E > E 2 however, the 
rate of type 3b regular orbits drops suddenly and for £ > £3 it vanishes. In the same vein, 
the percentage of escaping orbits is about 60% for £ = £2, while for larger values of energy 
it decreases reaching about one third of its initial value (20%) for £ = —1. The percentage 
of orbits that crash into the oblate primary on the other hand, grows drastically for £ > £2 
and at high values of energy (£ > £4) it seems to saturate to around 80%. Once more we 
observe, that the variation on the value of the total orbital energy has practically no influence 
on the rates of all the other types of orbits which remain unperturbed and at extremely low 
values (less than 5%). We seen in Fig. 12 that the (j) > 0 part of the configuration space 
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Fig. 12 The orbital structure of the physical (jr,y) plane in a corotating frame of reference is given using orbit 
type diagrams (OTDs) for three energy levels and for both parts (j) <0 (left column) and (j> >0 (right column) 
of the surface of section r = 0, when Aj = 0.1. (Top row): E — —1.90; (middle row): E — —1.70; (bottom 
row): E — —1.40. The vertical black dashed lines denote the centers of the two primaries, wile the vertical 
purple dashed lines indicate the position of the Lagrangian points L 2 and L 3 . The color bar contains the color 
code which relates the types of orbits presented in Fig.[^with different colors. 
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Fig. 13 Distribution of the escape and collisional time of the orbits on the physical (;c, y) space when A i =0.1 
for the energy levels of Fig. [T^ (Top row): E — —1.90; (middle row): E — —1.70; (bottom row): E — —1.40. 
The darker the color, the larger the escape/crash time. Initial conditions of bounded regular orbits are shown 
in white. 
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Fig. 14 Evolution of the percentages of all types of orbits on the physical (jr,}’) plane, for Ai = 0.1 when 
varying the energy parameter E for (a-left): the 0 < 0 part and (b-right): the 0 > 0 part of the configuration 
space. 


is dominated by initial conditions of orbits that escape from the system. Fig. [T4} t verifies 
this claim since the corresponding percentage remains high (more than 80%) throughout. 
Moreover one can observe that in this case crash on the oblate primary body becomes less 
and less unlike with increasing energy due to the fact that the rate of this family decreases 
as we proceed to higher energy levels and for £ = —1 they occupy only about 5% of the 
physical space. 

Before closing this section we would like to emphasize that the OTDs given in Figs. 
EE and [T2| have both fractal and non-fractal (smooth) boundary regions which separate 
the escape basins and the collisional basins. Such fractal basin boundaries is a common 
phenomenon in leaking Flamiltonian systems (e.g.,|Bleher et al. (|]^9^|l; JdeMoura^ Leteher] 


l |1999[ l; |de Moura & Grebogi| ( |2002[ l; [Schneider et al, (2002 2003 I; |Tuval et~ ( 2004| l). 


In the PCRTBP system the leakages are defined by both escape and crash conditions thus 
resulting in three exit modes. However, due to the high complexity of the basin boundaries, 
it is very difficult, or even impossible, to predict in these regions whether the test body (e.g., 
a satellite, asteroid, planet etc) collides with a primary body or escapes from the dynamical 
system. 


5 An overview analysis 

The color-coded OTDs in both parts of the physical {x,y) space provide sufficient informa¬ 
tion on the phase space mixing however, for only a fixed value of the energy integral and 
also for orbits that traverse the surface of section either directly (progradely) or retrogradely. 
Henon back in the late 60s ( |Henon[fl969[ l, introduced a new type of plane which can provide 
information not only about stability and chaotic regions but also about areas of bounded and 
unbounded motion using the section y = i = 0, y>0 (see also |Barrio et al.|f2008[ l). In other 
words, all the orbits of the test particles are launched from the x-axis with x = xq, parallel to 
the y-axis (y = 0). Consequently, in contrast to the previously discussed types of planes, only 
orbits with pericenters on the x-axis are included and therefore, the value of the energy E can 
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Fig. 15 (left column): Orbital structure of the {x,E) plane when (Top row): Ai — 0.001; (middle row): Aj = 
0 . 01 ; (bottom row): Ai = 0 . 1 . (right column): the distribution of the con'esponding escape/collisional times 
of the orbits. The vertical black dashed lines denote the centers of the two primaries, wile the vertical purple 
dashed lines indicate the position of the Lagrangian points L 2 and L 3 . The color bar contains the color code 
which relates the types of orbits presented in Fig.[^with different colors. 


be used as an ordinate. In this way, we can monitor how the energy influences the overall or¬ 
bital structure of our dynamical system using a continuous spectrum of energy values rather 
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Fig. 16 (left column): Orbital structure of the (x,Ai) plane when (Top row): E — —2.0; (middle row): E — 
— 1.5; (bottom row): E — —l.O. (right column): the distribution of the corresponding escape/collisional times 
of the orbits. The vertical black dashed lines denote the centers of the two primaries, wile the vertical purple 
dashed lines indicate the position of the Lagrangian points L 2 and L 3 . The color bar contains the color code 
which relates the types of orbits presented in Fig.[^with different colors. 


than few discrete energy levels. In the left column of Fig.[T^we present the orbital structure 
of the (v,£’) plane for three values of the oblateness coefficient when E E [—3,1], while in 
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the right column of the same figure the distribution of the corresponding escape/collision 
times of orbits is depicted. 

We observe the presence of several types of regular orbits around the two primary bod¬ 
ies. Being more precise, on both sides of the primaries we identify stability islands corre¬ 
sponding to both direct (counterclockwise) and retrograde (clockwise) quasi-periodic orbits. 
It is seen that a large portion of the exterior region, that is for x < x(L 2 ) and x > x{L-i), a 
large portion of the (jc, E ) plane is covered by initial conditions of escaping orbits however, 
at the left-hand side of the same plane two stability islands of type 3b regular orbits are 
observed. Additional numerical calculations reveal that for much lower values of x (jc<5) 
these two stability islands are joined and form a crescent-like shape. Furthermore, orbits 
with initial conditions very close to vertical line x = Ci, or in other words close to the center 
of the oblate primary 1 collide almost immediately with it, while their portion (thickness of 
the line) increases for larger values of the oblateness coefficient Ai. We also see that crash 
basins to primary 1 leak outside the interior region, mainly outside Li, and create compli¬ 
cated spiral shapes in the exterior region. On the other hand, the thin red bands represent 
initial conditions of orbits that collide with primary body 2. It should be pointed out that in 
the blow-ups of the diagram several additional very small islands of stability have been iden- 
tifie(0 We may say that the stability islands around primary body 2 are almost unperturbed 
by the shifting on the value of the oblateness. In contrast, we see that for A i =0.1 there is no 
indication of regular motion around oblate primary 1, while the amount of orbits that crash 
into primary 1 significantly grows with increasing Ai. Another interesting phenomenon is 
the fact that as primary body 1 becomes more and more oblate the fractility of the (x,E) 
plane reduces and the boundaries between escaping and collisional motion appear to be¬ 
come smoother. The following TablefTlshows the percentages of the all types of orbits in the 
color-coded OTDs shown in Figs.|15fand[T^ We observe that type 3a and 4 regular orbits 
are completely absent. It should be pointed out that the high values of escaping and crash 1 
percentages are due to the extended scattering region and large oblateness, respectively. 

Table 1 Percentages of all types of orbits in the color-coded OTDs shown in Figs. |l5[ a-c) and |l6[ a-c). ((a) 
corresponds to top row of the figures, (b) corresponds to middle rows, while (c) corresponds to bottom rows). 


Fij 

mre 

Escaping 

Crash 1 

Crash 2 

Type la 

Type lb 

Type 2a Type 2b 

Type 3a Type 3b 

Type 4 

Other 

15 

i 

64.29 

18.33 

1.46 

0.73 

2.04 

3.15 

1.41 

0.00 

6.93 

0.00 

1.66 

T5 


59.14 

27.53 

1.12 

0.13 

0.55 

3.09 

1.38 

0.00 

5.81 

0.00 

1.25 

T5 


52.03 

38.73 

0.86 

0.00 

0.01 

2.77 

1.06 

0.00 

4.08 

0.00 

0.46 

TC 

i 

28.54 

19.52 

1.71 

0.63 

1.08 

7.54 

5.15 

0.00 

35.82 

0.00 

0.01 

TC 


54.08 

36.09 

1.77 

0.00 

1.54 

6.47 

0.00 

0.00 

0.04 

0.00 

0.01 

TC 


40.31 

51.24 

2.05 

0.00 

2.53 

3.72 

0.00 

0.00 

0.00 

0.00 

0.15 


In order to obtain a more complete view of the orbital structure of the system, we follow 
a similar numerical approach to that explained before varying now the value of the oblate¬ 
ness coefficient Ai for three energy levels which belong to the three main Hill’s regions. 
This allow us to construct again a two-dimensional (2D) plane in which the x coordinate of 
orbits is the abscissa, while the logarithmic value of the oblateness coefficient log[Q(Ai) is 
the ordinate. The orbital structure of the (x,Ai) plane when log[g(Ai) e [-3,-1] is shown 
in the left column of Fig. m while the distribution of the corresponding escape/collision 

^ An infinite number of regions of (stable) quasi-periodic (or small scale chaotic) motion is expected from 
classical chaos theory. 
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times of orbits is given in the the right column of the same figure. The black solid line is the 
limiting curve which distinguishes between regions of allowed and forbidden motion and is 
defined as 

f{y,Ai)=V{x = Q,y-Ai) = E, (7) 

while the vertical dashed magenta lines indicate the position of the Lagrangian points L 2 
and L 3 . Here it should be pointed out that as it was mentioned in Section]^ the position of 
the Lagrangian points L 2 and Lj is no longer fixed with variable Ai, since it is function of 
the oblateness coefficient. 

A very complicated orbital structure is reveled in the (x, Ai) plane from which however, 
we can deduce some interesting results such as: (i) for Z? = — 2 < £2 we see a very organized 
structure in which all the different domains, corresponding to different types of orbits, are 
well-defined thus the fractility of the plane is minimum, while for larger values of the energy 
the boundaries of the domains become highly fractal; (ii) orbits with initial conditions very 
close to the center of the oblate primary 1 collide with it practically immediately; (iii) the 
amount of collisional orbits to primary 1 grows with increasing energy and oblateness, while 
that of escaping orbits seems to reduces; (iv) the extent of the area on the (x,Ai) plane 
occupied by retrograde (clockwise rotating) orbits is reduced with increasing oblateness 
and is favoured only at relatively high energy levels (Z? > E^)\ (v) the stability island of 
retrograde motion around primary 2 is present only at low energies {E < E 2 ) and at low 
values of A\ {A\ < 0.01), while for larger values of energy or/and oblateness disappears; 
(vi) the stability island of prograde (direct) motion around primary 2 seems to be unaffected 
by the change on the value of the oblateness however, it slightly reduces as we proceed to 
higher energy levels. The phenomenon that stability islands can appear and disappear as a 
dynamical parameter is changed has also been reported in earlier paper (e.g., [Barrio et al.| 
( |2006| l; [de Assis & Terra| ( [2014| )). 


6 Discussion and conclusions 

The scope of this work was to shed some light to the properties of motion in the Copenhagen 
problem where one of the primaries is an oblate spheroid. We continued the work initiated in 
Paper I and II following similar numerical techniques therefore, this paper should be consid¬ 
ered as Part III. We managed to distinguish between bounded, escaping and collisional orbits 
and we also located the basins of escape/collison, finding correlations with the correspond¬ 
ing escape/crash times of the orbits. Our extensive and thorough numerical investigation 
strongly suggests, that the overall motion of a test body under the gravitational field of two 
primaries is a very complicated procedure and very dependent on the value of the energy 
integral. To our knowledge, this is the first detailed and systematic numerical analysis on the 
phase space mixing of bounded motion, escape and crash in the Copenhagen problem with 
oblateness and this is exactly the novelty and the contribution of the current work. 

We defined for several values of the total orbital energy E, dense uniform grids of 1024 x 
1024 initial conditions regularly distributed on both parts {(j) < 0 and (j) > 0) the physical 
(x,y) plane inside the area allowed by the value of the energy. All orbits were launched 
with initial conditions inside the scattering region, which in our case was a square grid 
with —2 < x,y < 2. For the numerical Integration of the orbits in each type of grid, we 
needed about between 7 hours and 8 days of CPU time on a Pentium Dual-Core 2.2 GHz 
PC, depending on the escape and collisional rates of orbits in each case. For each initial 
condition, the maximum time of the numerical integration was set to be equal to 10 "* time 
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units however, when a particle escaped or collided with one of the primaries the numerical 
integration was effectively ended and proceeded to the next available initial condition. 

The present article provides quantitative information regarding the escape and crash 
dynamics in the Copenhagen problem with oblateness. The main numerical results of our 
research can be summarized as follows: 

1. The value of the oblateness coefficient of the primary body 1 was found to greatly in¬ 
fluence the stability regions around the same primary body. In particular, as the primary 
becomes more and more oblate the size of the stability islands corresponding to regular 
motion around only primary 1 is reduced and at relatively high values of the oblateness 
there is no indication of regular motion around the oblate primary. 

2. We found that orbits with initial conditions in the vicinity of the oblate primary body 1 
crash almost immediately on it, while on the other hand orbits initiated around spherical 
primary 2 have larger non zero crash times. 

3. In both parts of the configuration space ((^ < 0 and (j) > 0) the crash basins of oblate 
primary 1 were found to be well defined broad and extended regions, while the initial 
conditions leading to crash with spherical primary 2 form thin filaments and spiral bands. 

fudging by the detailed and novel outcomes we may say that our task has been success¬ 
fully completed. We hope that the present numerical analysis and the corresponding results 
to be useful in the field of escape dynamics in the Copenhagen problem with oblateness. 
The outcomes as well as the conclusions of the present research are considered, as an initial 
effort and also as a promising step in the task of understanding the escape mechanism of 
orbits in this interesting version of the classical three-body problem. Taking into account 
that our results are encouraging, it is in our future plans to properly modify our dynamical 
model in order to expand our investigation into three dimensions and explore the entire six¬ 
dimensional phase thus revealing the influence of the oblateness coefficient on the orbital 
structure. Moreover, it would be interesting to apply our numerical methods in some inter¬ 
esting cases of the PCRTBP such the Earth-Moon and Saturn-Titan systems where however 
the values of oblateness are much smaller than those considered in the present paper. 
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